% !TEX root = FDS_Validation_Guide.tex

\chapter{Quantifying Model Uncertainty}
\label{Error_Chapter}

This chapter describes a method to estimate the model uncertainty using comparisons of model predictions with experimental measurements whose uncertainty has been quantified. The method is ideal for complex numerical models like FDS for which a systematic analysis of sub-components is impractical, but for which there exists a relatively large amount of experimental data with which to evaluate the accuracy of the model predictions. If the uncertainty in the experiments can be quantified, the uncertainty in the model can then be expressed in the form of a normal distribution whose mean and standard deviation are estimated from the relative difference between the predicted and measured values.

This method only addresses model uncertainty. It does not account for the uncertainty associated with the model input parameters. How the {\em parameter uncertainty} is treated depends greatly on the type of application. Regardless of how the parameter uncertainty is calculated, the model uncertainty needs to be addressed independently. In fact, it is incumbent on the model developers to express the uncertainty of the model in as simple a form as possible to enable the end user to assess the impact of parameter uncertainty and then combine the two forms of uncertainty into a final result.


\section{Introduction}

The most effective way of introducing the subject of uncertainty in fire modeling is by way of an example. Suppose that a fire model is used to estimate the likelihood that an electrical control cable could be damaged by a fire. It is assumed that the cable loses functionality when its surface temperature reaches 200~$^\circ$C, and the model predicts that the cable temperature could reach as high as 175~$^\circ$C. Does this mean that there is no chance of damage? The answer is no, because the input parameters, like the heat release rate of the fire, and the model assumptions, like the way the cables are modeled, are uncertain. The combination of the two -- the {\em parameter uncertainty} and the {\em model uncertainty} -- leave open the possibility that the cable temperature could exceed 200~$^\circ$C.

This chapter addresses {\em model uncertainty} only and suggests a simple method to quantify it. While parameter uncertainty is certainly an
important consideration in fire modeling, its treatment varies considerably depending on the particular application. For example, in what is often
referred to as a ``bounding analysis,'' the model input parameters are chosen so as to simulate a ``maximum credible'' or ``worst case'' fire. In
other cases, mean values of the input parameters constitute a ``50th percentile'' design scenario. Sometimes entire statistical distributions, rather
than individual values, of the input parameters are ``propagated'' through the model in a variety of ways, leading to a statistical distribution of
the model output. Notarianni and Parry survey these techniques in the SFPE Handbook~\cite{Notarianni:SFPE}. Regardless of the method that is chosen
for assessing the impact of the input parameters on the model prediction, there needs to be a way of quantifying the uncertainty of the model itself.
In other words, how good is the prediction for a given set of input parameters?

The issue of model uncertainty has been around as long as the models themselves. The scenario above, for example, was considered by Siu and Apostolakis in the early 1980s~\cite{Siu:RE1982} as part of their development of risk models for nuclear power plants. The fire models at the time were relatively simple. In fact, many were engineering correlations in the form of simple formulae. This made the methods for quantifying their uncertainty reasonably tractable because each formula consisted of only a handful of physical parameters. Over the past thirty years, however, both fire modeling and the corresponding methods of uncertainty analysis have become far more complex. The current generation of computational fluid dynamics (CFD) based fire models require such a large number of physical and numerical parameters that it is considered too cumbersome to estimate model uncertainty by systematically assessing their combined effect on the final prediction. The more practical approach is to compare model predictions with actual fire experiments in a validation study, the conclusions of which typically come in the form of statements like: ``The model generally over-predicts the measured temperatures by about 10~\%,'' or ``The model predictions are within about 20~\% of the measured heat fluxes.'' This information is helpful, at the very least to demonstrate that the model is appropriate for the given application. However, even the statement that the model over-predicts measured temperatures by 10~\% is useful not only to gain acceptance of the model but also to provide a better sense of the model's accuracy, and a greater level of assurance in answering the question posed above. Knowing that the model not only predicted a temperature of 175~$^\circ$C, but also that the model tends to over-predict temperatures by a certain amount, increases the confidence that the postulated fire would not cause the cable to fail. The probability of cable failure could be quantified further via a statistical distribution like the one shown in Fig.~\ref{bell_curve}. The area indicated by the shaded region is the probability that the temperature will exceed 200~$^\circ$C, even though the model has predicted a peak temperature of only 175~$^\circ$C.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=5.in]{FIGURES/bell_curve}
\end{center}
\caption[Demonstration of model uncertainty]{Plot showing a possible way of expressing the uncertainty of the model prediction.}
\label{bell_curve}
\end{figure}

This chapter describes a method for expressing {\em model uncertainty} by way of a distribution like the one shown in Fig.~\ref{bell_curve}. The procedure is not a dramatic departure from the current practice of fire model validation in that it relies entirely on comparisons of model predictions and experimental measurements. The advantage of the approach is that it does not demand advanced knowledge of statistics or details of the numerical model. The parameters of the distribution shown in Fig.~\ref{bell_curve}, namely the mean and standard deviation, are not generated by the model user. Rather, they are reported as the results of the validation study. The calculation of the probability of exceeding some critical threshold (i.e., the area under the curve) is a simple table look-up or function call in data analysis software like Microsoft Excel\textregistered.



\section{Sources of Model Uncertainty}

A deterministic fire model is based on fundamental conservation laws of mass, momentum and energy,
applied either to entire compartments or smaller control
volumes that make up the compartments. A CFD model may use millions of control volumes to compute the
solution of the Navier-Stokes equations.
However, it does not actually solve the Navier-Stokes equations, but rather an approximate form of these equations.
The approximation involves simplifying
physical assumptions, like the various techniques for treating subgrid-scale turbulence.
One critical approximation is the discretization of the governing equations. For example,
the partial derivative of the density, $\rho$,
with respect to the spatial coordinate, $x$, can be written in approximate form as:
\be \frac{\partial \rho}{\partial x} = \frac{\rho_{i+1} - \rho_{i-1}}{2 \, \dx} + \mathcal{O}(\dx^2) \ee
where $\rho_i$ is the average value of the density in the $i$th grid cell and $\dx$ is the spacing between cells.
The second term on the right represents all of the terms of order $\dx^2$ and higher in the Taylor
series expansion and are known collectively as the
{\em discretization error}. These extra terms are simply dropped from
the equation set, the argument being that they become smaller and smaller with decreasing grid cell size, $\dx$.
The effect of these neglected terms is captured, to
some extent, by the subgrid-scale turbulence model, but that is yet another approximation of the true physics.
What effect do these approximations have on
the predicted results? It is very difficult to determine based on an analysis of the discretized equations.
One possibility for estimating
the magnitude of the discretization error is to perform a detailed
convergence analysis, but this still does not answer a
question like, ``What is the uncertainty of the model prediction of the gas
temperature at a particular location in the room at a particular point in time?''

To make matters worse, there are literally dozens of subroutines that make up a CFD-based fire model,
from its transport equations, radiation solver, solid phase heat transfer routines, pyrolysis model,
empirical mass, momentum and energy transfer routines at the wall, and so on.
It has been suggested by some that
a means to quantify the model uncertainty is to combine the uncertainties of all the model
components.
However, such an exercise is very difficult, especially for a CFD model,
for a number of reasons. First, fire involves
a complicated interaction of gas and solid phase phenomena that are closely coupled.
Second, grid sensitivity in a CFD model or the error associated with
a two-layer assumption in a zone model are dependent on the particular fire scenario.
Third, fire is an inherently transient phenomenon in which relatively small
changes in events, like a door opening or sprinkler actuation, can lead to significant changes in outcome.

Rather than attempt to decompose the model into its constituent parts and assess the uncertainty of
each, the strategy adopted here is to compare model predictions to as many
experiments as possible. This has been the traditional approach for quantifying model uncertainty in fire
protection engineering because of the relative abundance of test data. Consider, for example, the
plot shown in Fig.~\ref{scatterplot}. This is the typical outcome of a validation study, where in this case a series of
heat flux measurements are compared with model predictions.
The diagonal line indicates where the prediction and measurement agree.
But because there is uncertainty associated with each, it cannot be said that the model is perfect if its predictions
agree exactly with measurements.
There needs to be a way of quantifying the uncertainties of each before any conclusions can be drawn.
Such an exercise would result in the uncertainty
bars\footnote{The data in Fig.~\ref{scatterplot} was extracted from Ref.~\cite{NUREG_1824_Sup_1}.
The uncertainty bars are for demonstration only.}
shown in the figure. The
horizontal bar associated with each point represents the uncertainty in the measurement itself.
For example, the heat flux gauge is subject to uncertainty due to its design and fabrication.
Because the horizontal bar represents the experimental uncertainty, it is assumed that the vertical
bar represents the model uncertainty. This is only partially true. In fact the vertical bar represents the total
uncertainty of the prediction, which is a combination of the {\em model} and {\em parameter} uncertainties. The physical
input parameters, like the heat release rate and material properties, are based on measurements that are reported
in the documentation of the experiment.
The total {\em experimental uncertainty} is represented by all of the horizontal bar and part of the vertical.
If the {\em experimental uncertainty} can be quantified, then the {\em model uncertainty} can be obtained as a result.


\begin{figure}[ht]
\begin{center}
\includegraphics[height=3.in]{FIGURES/scatterplot}
\end{center}
\caption[Sample scatter plot]{Example of a typical scatter plot of model predictions and experimental measurements.}
\label{scatterplot}
\end{figure}


\section{Experimental Uncertainty}

The difference between a model prediction and an experimental measurement is a combination of three components: (1) uncertainty in the measurement of the predicted quantity, (2) uncertainty in the model input parameters, and (3) uncertainty in the model physics and numerics. The first two components are related to uncertainty in the measured input and output quantities. For example, consider the hot gas layer (HGL) temperature. First, the thermocouple measurements used to calculate the HGL temperature have uncertainty. Second, the measurement of the heat release rate of the fire has uncertainty, and this uncertainty affects the predicted HGL temperature. Third, the model itself, including its physical assumptions and numerical approximations, has uncertainty. The objective of the validation study is to quantify this third component, the model uncertainty. To do this, the first two components of uncertainty related to the experimental measurements must be quantified. The combination of these two, the experimental uncertainty, is the objective of this section.

For many of the experiments considered in this guide, the uncertainty of the measurements was not documented. Instead, estimates of measurement uncertainty are made based on those few experiments that do include uncertainty estimates, and this information is supplemented by engineering judgment. In the following two subsections, each component of the experimental uncertainty is considered separately. First, the uncertainty in the measurement of the predicted quantity of interest, like the surface temperature of the compartment. Second, the uncertainties of the most important input parameters are propagated through simple models to quantify their effect on the predicted quantity. Then, the uncertainties are combined via simple quadrature to estimate the total experimental uncertainty.

Note that in this guide, all uncertainties are expressed in relative form, as a percentage. The uncertainty of a measurement is a combination of the systematic uncertainty associated with the various underlying measurements and assumptions; and the random uncertainty associated with the conduct of the experiment. Following the recommended guidelines for evaluating and expressing the uncertainty of measurements~\cite{Taylor&Kuyatt:1994}, the systematic and random uncertainty values are combined via quadrature resulting in a combined relative standard uncertainty.

\subsection{Uncertainty of Common Fire Measurements}

Because most of the experiments described in this guide have little or no information about the uncertainty of the measurements, much of this section is based on the uncertainty analysis contained in the test report of the NIST/NRC Experiments\footnote{Note that the uncertainties in Ref.~\cite{Hamins:SP1013-1} are reported in the form of 95~\% confidence intervals, or ``2-sigma''. This is referred to as the expanded uncertainty with a coverage factor of 2. To avoid confusion, in this report the uncertainty of all measurements and model predictions shall be reported in terms of a relative standard uncertainty; that is, the estimated standard deviation of the measured or predicted quantity.}. The types of measurements described in this report are the ones most commonly used in large scale fire experiments. They include thermocouples for gas and surface temperature measurements, heat flux gauges, smoke and gas analyzers, and pressure sensors.

\subsubsection{Thermocouples}

Thermocouples are used to measure both gas and surface temperatures. They come in a variety of sizes and are constructed of different types of metals. Some are shielded or aspirated to limit the influence of thermal radiation from remote sources. In Ref.~\cite{Hamins:SP1013-1}, Hamins et al. estimate the uncertainty of the various thermocouple measurements. Estimates of the combined relative standard uncertainty fall in a range between 2.5~\% and 7.5~\%. Because it is not possible to analyze the thousands of thermocouple measurements made in the experiments, the relative standard uncertainty applied to all thermocouple measurements is 5~\%.

\subsubsection{Heat Flux Gauges}

For the NIST/NRC experiments, four types of heat flux gauges were used, some of which measured the total heat flux, and some of which measured only the radiation heat flux. The uncertainty associated with a heat flux measurement depends on many factors, including gauge characteristics, the calibration conditions and accuracy, as well as the incident flux modes (convective, radiative, conductive). Typically, the reported relative standard uncertainty of heat flux gauges varies from about 2.5~\% to 5~\%, with the measurement uncertainty dominated by uncertainty in the calibration and repeatability of the measurement. Repeatability of the various heat flux measurements in the NIST/NRC experiments was determined by examining measurements by the same instruments for different pairs of repeat tests. The difference between the measurements was about 3.5~\%, on average, for both the radiative flux measurements and the total flux measurements. For all of the experiments described in this guide, a combined relative standard uncertainty of 5~\% is suggested based on the measurement repeatability and calibration uncertainties reported for the NIST/NRC experiments.

\subsubsection{Gas Analyzers}

Gas concentrations were measured in two sets of experiments conducted at NIST, the NIST/NRC and the WTC experiments. The volume fractions of the combustion products, carbon monoxide (CO) and carbon dioxide (CO$_2$), were measured using gas sampling in conjunction with non-dispersive infrared analyzers, while the oxygen (O$_2$) volume fraction was typically measured using a paramagnetic analyzer. Gases were extracted through stainless steel or other types of lines and were pumped from the compartment and passed through the analyzers. For several reasons, water in the sample was typically filtered, so the reported results are denoted as ``dry''. Analyzers were calibrated through the use of standard gas mixtures, with low relative uncertainties. Problems with the technique may involve instrument drift, analyzer response, incomplete and partial drying of sample gases, or (in the case when drying is not used) undetermined amounts of water vapor in the oxygen cell, which result in inaccurate readings.

For the NIST/NRC experiments, the species were measured in both the upper and lower layers. The relative standard uncertainty in the measured values was about 1.5~\% for both the O$_2$ depletion and the CO$_2$ measurements. The largest contributors were the uncertainty in the composition of the calibration gas and the possibility of an undetermined amount of water vapor in the sample. The difference between the repeat measurements was about 1~\%, on average, for both the O$_2$ depletion and the CO$_2$ measurements. Therefore, the combined relative standard uncertainty is estimated to be 2~\% for measurements of both the O$_2$ decrease and the CO$_2$ increase.

\subsubsection{Smoke Light Extinction Calculations}

The smoke concentration was measured in the NIST/NRC experiments using laser transmission at 632.8~nm. The reported mass concentration of smoke, $m_s'''$, was computed using the following expression:
\be
   m_s'''=\frac{\ln(I_0/I)}{\phi_{\rm s} L}
\ee
where $L$ is the path length, $I$ and $I_0$ are the laser signal and reference signal, respectively, and $\phi_{\rm s}$ is the specific extinction coefficient, which has a nearly universal value of 8.7~m$^2$/g~$\pm$~5~\% for hydrocarbons~\cite{Mulholland:F+M}. The systematic relative standard uncertainty of this measurement was reported to be 9~\% in Ref.~\cite{Hamins:SP1013-1}, with the dominant contribution to the uncertainty coming from drift in the laser measurement. Repeatability of the smoke measurement was investigated for the NIST/NRC experiments. The mean difference between replicate measurements was about 11~\%. Therefore, a combined relative standard uncertainty of 14~\% is suggested.

\subsubsection{Pressure Gauges}

The uncertainty in pressure measurements is typically small, but depends on the sensor type and calibration. In the NIST/NRC experiments, the differential pressure gauge used was temperature compensated, highly linear, and very stable. The estimated relative standard uncertainty is 0.5~\%.

\subsubsection{Bi-Directional Probes}

Gas velocity is typically measured in fire experiments using bi-directional probes. These devices work like pitot tubes but have much larger openings. Bryant~\cite{Bryant:FSJ2009} estimates that the standard relative uncertainty of this measurement, assuming that the probe is aligned well with the flow, is approximately 7~\%.

\subsubsection{Oxygen Consumption Calorimeters}

For all of the experiments described in this guide, the heat release rate (HRR) is determined either via oxygen consumption calorimetry or via the mass loss rate multiplied by the fuel heat of combustion. The accuracy of each method varies roughly between 2.5~\%, where the fire is small and the fuel stoichiometry is well understood, and 13~\%, where the fire is large or the smoke is not completely captured or the fuel stoichiometry is not well understood. In Ref.~\cite{Hamins:SP1013-1}, the relative standard uncertainty of a 2~MW heptane spray fire is estimated to be 7.5~\%. It is assumed that the uncertainty of the HRR for the other experiments is comparable.

\subsubsection{Device Activation or Failure Time}

Fire models are often used to predict the time to activation of devices like sprinklers and smoke detectors, and time to failure of critical equipment like electrical cables. Measuring activation or failure times in experiments is fairly precise, and, thus, the uncertainty of such measurements is essentially zero. Almost all of the uncertainty associated with these times is in the measurement or characterization of the mechanism of activation or failure. For example, the activation of a sprinkler is based on its measured RTI (Response Time Index) and activation temperature. Estimates of the uncertainty of these parameters are discussed in the next section.

\subsection{Propagation of Input Parameter Uncertainty}

Empirical correlations for basic fire phenomena provide a convenient way to assess the propagation of the uncertainty of the model input parameters. The more complex fire models may require dozens of physical and numerical input parameters for a given fire scenario.  However, only a few of these parameters, when varied over their plausible range of values, significantly impact the results. For example, the thermal conductivity of the compartment walls does not significantly affect a predicted cable surface temperature, but the HRR of the fire does. The relatively simple empirical models identify the key parameters that impact the predicted quantity, and they provide the means to quantify the functional relationship between model inputs and outputs.

\subsubsection{Gas and Surface Temperatures}

According to the McCaffrey, Quintiere, and Harkleroad (MQH) correlation, the HGL temperature rise, $T-T_0$, in compartment fire is proportional to the HRR, $\dot{Q}$, raised to the two-thirds power:
\be
   T-T_0 = C \, \dot{Q}^{\frac{2}{3}}
\ee
The constant, $C$, involves a number of geometric and thermo-physical parameters that are unique to the given fire scenario. By way of differentials, this empirical relationship can be expressed in the form:
\be
   \frac{\Delta T}{T-T_0} \approx \frac{2}{3} \, \frac{\Delta \dot{Q}}{\dot{Q}} \label{MQH_diff}
\ee
This is a simple formula with which one can readily estimate the relative change in the temperature rise due to the relative change in the HRR. The uncertainty in the HRR of the validation experiments was estimated to be 7.5~\%. Equation~(\ref{MQH_diff}) indicates that a 7.5~\% increase in the HRR should lead to a 5~\% increase in the HGL temperature.

\subsubsection{HGL Depth}

Most of the experiments for which the HGL depth was predicted had at least one open door or window that effectively determined the steady-state HGL depth. Unlike all of the other predicted quantities, the HGL depth is relatively insensitive to the fire's HRR. It is largely determined by the height of the opening, and for this reason there is essentially no uncertainty associated with the model inputs that affect the layer depth.

\subsubsection{Gas and Smoke Concentration}

Most fire models assume that combustion product gases and soot, once beyond the flaming region of the fire, are passively transported throughout the compartment. The major products of combustion, like CO$_2$ and water vapor, plus the major reactant, O$_2$, are generated, or consumed, in direct proportion to the burning rate of the fuel, which is directly proportional to the HRR. The mass fraction of any species in the HGL is directly proportional to the product of its yield and the HRR.

For many of the experiments described in this guide, the yields of the major product gases like O$_2$ and CO$_2$ from pure fuels like methane gas and heptane liquid are known from the basic stoichiometry to a high level of accuracy. Thus, the relative uncertainty in the concentration of major products gases is the same as that of the HRR, 7.5~\%. The uncertainty in the smoke concentration, however, is a combination of the uncertainty of the HRR and the soot yield. The relative standard uncertainty of the soot yield of heptane reported in Ref.~\cite{Hamins:SP1013-1} is 11~\%. The uncertainties for HRR and soot yield are combined via quadrature and the resulting expanded relative uncertainty is 13~\%.

\subsubsection{Pressure}

In a closed, ventilated compartment, the average pressure, $p$, is governed by the following equation:
\be
   \frac{\rm{d} p}{\rm{d} t} = \frac{\gamma-1}{V} \left( \dot{Q} - \dot{Q}_{\rm loss} \right) + \frac{\gamma p}{V} \left( \dot{V}-\dot{V}_{\rm leak} \right) \label{simple_pressure}
\ee
where $\gamma$ is the ratio of specific heats (about 1.4), $V$ is the compartment volume, $\dot{Q}$ is the HRR, $\dot{Q}_{\rm loss}$ is the sum of all heat losses to the walls, $\dot{V}$ is the net ventilation rate into the compartment, and $\dot{V}_{\rm leak}$ is the leakage rate out of the compartment. The leakage rate is a function of the compartment over-pressure:
\be
   \dot{V}_{\rm leak} = A_{\rm leak} \sqrt{ \frac{2 (p-p_\infty)}{\rho_\infty} }
\ee
The maximum compartment pressure is achieved when the pressure rise term in Eq.~(\ref{simple_pressure}) is set to zero. Rearranging terms yields an estimate for the maximum pressure:
\be
   (p-p_\infty)_{\max} \approx \frac{\rho_\infty}{2} \left( \frac{ (\gamma-1)  \left( \dot{Q} - \dot{Q}_{\rm loss} \right) + \gamma p_\infty \dot{V} }{\gamma p_\infty A_{\rm leak} } \right)^2 \label{pmax}
\ee
The test report for the NIST/NRC experiments contains estimates of the uncertainty in the HRR, ventilation rate and leakage area. To calculate the uncertainty in the maximum pressure rise resulting from the uncertainty in these three parameters, the pressure rise estimate in Eq.~(\ref{pmax}) was calculated using 1000 randomly selected sets of values of the HRR, ventilation rate, and leakage area. These parameters were assumed to be randomly distributed with mean values of 1000~kW, 1~m$^3$/s, and 0.06~m$^2$ and relative standard uncertainties of 75~kW, 0.1~m$^3$/s, and 0.0021~m$^2$. The mean values of these parameters were typical of the NIST/NRC experiments, and the uncertainties were reported in the test report. The resulting relative standard uncertainty in the pressure of a closed compartment due to the uncertainty in the HRR, ventilation rate, and leakage area is 21~\%.

For an open compartment, in which the ventilation rate and leakage area have much less influence, the relative standard uncertainty in the pressure is twice that of the HRR, 15~\%.

\subsubsection{Velocity}

Fire-induced velocities, as in a plume or ceiling jet, are roughtly proportional to the HRR to the 1/3 power~\cite{SFPE:Alpert}. Given that the relative standard uncertainty in the HRR is 7.5~\%, the uncertainty in gas velocity due to the propagated effect of the uncertainty in the HRR is 2.5~\%.

\subsubsection{Heat Flux}

The heat flux to a target or wall is a combination of direct thermal radiation from the fire and convective and thermal radiation from the HGL. If the heat flux is predominantly due to the thermal radiation of the fire, it can be approximated using the point source radiation model:
\be
   \dot{q}'' = \frac{ \chi_{\rm r} \dot{Q} }{4 \pi r^2}
\ee
where $\chi_{\rm r}$ is the radiative fraction, $\dot{Q}$ is the HRR, and $r$ is the distance from the fire. The relative standard uncertainty of the heat flux is a combination of the uncertainty in the radiative fraction and the HRR:
\be
   \frac{ \delta \dot{q}''}{\dot{q}''} \approx \frac{\delta \dot{Q}}{\dot{Q}} + \frac{\delta \chi_{\rm r}}{\chi_{\rm r}}
\ee
Reference~\cite{Hamins:SP1013-1} estimates the relative standard uncertainty of the radiative fraction of a heptane pool fire to be 8~\%. Combined with the 7.5~\% uncertainty in the HRR (via quadrature) yields a 11~\% relative standard uncertainty in the heat flux directly from a fire.

The heat flux to a cold surface due to the exposure to hot gases and not necessarily the fire itself is the sum of radiative and convective components:
\be
   \dot{q}'' = \epsilon \sigma \left( T_{\rm gas}^4 - T_\infty^4 \right) + h \left( T_{\rm gas} - T_\infty \right)
\ee
where $\epsilon$ is the surface emissivity, $\sigma$ is the Stefan-Boltzmann constant, $T_{\rm gas}$  is the gas temperature, $T_\infty$ is the ambient temperature, and $h$ is the convective heat transfer coefficient. From the discussion above, the relative standard uncertainty in the gas temperature rise above ambient is 5~\% resulting from an estimated uncertainty in the HRR of 7.5~\%. There is also uncertainty in the convective heat transfer coefficient, but this is attributed to the model, not the experimental measurements. Thus, the uncertainty in the heat flux is largely a function of the uncertainty in the gas temperature which is largely a function of the HRR. As was done for the pressure, 1000 randomly selected values of gas temperature with a mean of 300~$^\circ$C above ambient and an relative uncertainty of 5~\% resulted in a corresponding uncertainty of 9~\% in the heat flux.

In actual compartment fires, the heat flux to surfaces is a combination of direct thermal radiation from the fire and indirect radiation and convection from the hot gases. Given that the calculation of the former incurs a 11~\% relative standard uncertainty and the latter 9~\%, to simplify the analyses, a value of 10~\% is used for all heat flux predictions.

\subsubsection{Sprinkler Activation Time}
\label{uncertainty_sprinkler_acts}

The uncertainty in the reported sprinkler activation times is due mainly to uncertainties in the measured HRR, RTI (Response Time Index), and activation temperature. There is a negligible uncertainty in the measured activation time itself, which is typically determined with a pressure transducer. To determine the effect of the uncertainties in the HRR, RTI and activation temperature, consider the ordinary differential equation governing the temperature, $T_{\rm link}$, of a conventional glass bulb of fusible link sprinkler:
\be
   \frac{\d T_{\rm link}}{\d t} = \frac{\sqrt{u}}{\rm RTI} \left( T_{\rm gas} - T_{\rm link} \right) \label{RTIeq}
\ee
Here, $u$ and $T_{\rm gas}$  are the velocity and the temperature of the ceiling jet, respectively. According to Alpert's ceiling jet correlation~\cite{SFPE:Alpert}, the ceiling jet temperature and velocity are proportional to the HRR raised to the power of 2/3 and 1/3, respectively. Given the relative standard uncertainty in the HRR of 7.5~\%, the uncertainty in the ceiling jet temperature and velocity are, thus, 5~\% and 2.5~\%, respectively. As for the RTI and activation temperature, these values are measured experimentally and the uncertainties differ depending on the test procedure. Vettori~\cite{Vettori:1} reports that the RTI of the sprinklers used in his experiments is 56~(m$\cdot$s)$^{1/2}$ with a relative standard uncertainty of 11~\%, and that the activation temperature is 68~$^\circ$C~$\pm$~2.4~$^\circ$C. This latter uncertainty estimate is assumed to represent one standard deviation. Assuming an ambient temperature of approximately 20~$^\circ$C, the relative standard uncertainty in the activation temperature is assumed to be 5~\%.

To determine how the uncertainties in the measured parameters affect the sprinkler activation time, Eq.~(\ref{RTIeq}) was integrated 1000 times using random selections of the ceiling jet temperature and velocity, RTI, and activation temperature. The mean ceiling jet temperature was increased linearly at rates varying from 0.5~$^\circ$C/s to 2~$^\circ$C/s, consistent with the variety of growth rates measured by Vettori. The mean ceiling jet velocity was assumed to be 1~m/s. This procedure yielded a relative standard uncertainty in the sprinkler activation time of 6~\%.

The activation times recorded by Vettori include two or three replicates for each configuration. The standard deviation of the 45 measured activation times, normalized by the mean of each set of replicates, was 6~\%, consistent with the result obtained above.

\subsubsection{Number of Activated Sprinklers}

Alpert's ceiling jet correlation~\cite{SFPE:Alpert} predicts the temperature rise, $T-T_\infty$ ($^\circ$C), as a function of the HRR, $\dot{Q}$ (kW), and radial distance, $r$ (m), from the plume centerline:
\be
   T-T_\infty = 5.38 \; \frac{ \dot{Q}^{2/3} / H^{5/3} }{ (r/H)^{2/3} } \quad ; \quad r/H > 0.18
\ee
For a given ceiling height, $H$, the radial extent of the sprinkler activation temperature is directly proportional to $\dot{Q}$. The number of activated sprinklers is roughly proportional to the square of this radial distance, assuming the sprinklers are uniformly spaced on a rectangular grid. Thus, the uncertainty in the number of activated sprinklers due to the uncertainty in the HRR is 15~\%.


\subsubsection{Electrical Cable Failure Time}

The uncertainty in the reported cable failure times is due mainly to uncertainties in the measured exposing temperature, cable diameter, and jacket thickness. The uncertainty in the measured mass per unit length of the cable is assumed to be negligible. To determine the uncertainty in the cable failure time, the heat conduction equation in the THIEF model was solved numerically using 10,000 random selections of the exposing temperature, cable diameter, and jacket thickness. The cable diameter was varied from 16.25~mm to 16.35~mm, and the jacket thickness was varied from 1.45~mm to 1.55~mm. The uncertainty in the exposing temperature of the cylindrical heater was assumed to be 2.5~\%, the lower bound of the range of uncertainty estimates for thermocouple measurements. The mass per unit length of the cable was assumed to be 0.529~kg/m, and the ambient temperature was assumed to be 20~$^\circ$C. This procedure yielded an estimated relative standard uncertainty in the cable failure time of 12~\%.

\subsubsection{Smoke Detector Activation Time}

There is a single set of experiments with which to evaluate model predictions of smoke detector activation time, the NIST Home Smoke Alarm Experiments. The test report~\cite{Bukowski:1} does not include detailed information about the alarm mechanism within the various smoke detectors used in the experiments. Thus, from a modeling standpoint, these devices are ``black boxes'' and their activation can only be discerned from a variety of empirical techniques, the most popular of which is to assume that the smoke detector behaves like a sprinkler or heat detector whose activation is governed by Eq.~(\ref{RTIeq}) with a low activation temperature and RTI. Bukowski and Averill~\cite{Bukowski:2} suggest an activation temperature of 5~$^\circ$C to be typical of many residential smoke alarms. The propagated uncertainty of this estimate is difficult to determine because temperature rise is not particularly well-correlated with smoke concentration within the sensing chamber of the detector. Nevertheless, the relative standard deviation of the normalized activation times for the NIST Home Smoke Alarm Experiments is 34~\%. Without more detailed information about the activation criteria, the models cannot predict the activation times more accurately than this value.

\subsection{Summary of Experimental Uncertainty Estimates}

Table~\ref{Uncertainty} summarizes the estimated uncertainties of the major output quantities. The right-most column in the table represents the total experimental uncertainty, denoted as $\widetilde{\sigma}_E$, a combination of the uncertainty in the measurement of the output quantity itself, along with the propagated uncertainties of the key measured input quantities. This total experimental uncertainty is obtained by taking the square root of the sum of the squares of the measurement and propagation uncertainties that have been estimated in the previous two sections. It is assumed that the two forms of uncertainty are independent.

\begin{table}[ht]
\caption[Summary of uncertainty estimates]{Summary of uncertainty estimates. All values are expressed in the form of a standard relative uncertainty.}
\begin{center}
\begin{tabular}{|l|c|c|c|}
\hline
Output                        & Measurement     & Propagated Input  & Combined                                  \\
Quantity                      & Uncertainty     & Uncertainty       & Uncertainty, $\widetilde{\sigma}_E$       \\ \hline \hline
Gas and Solid Temperatures    & 0.05            & 0.05              & 0.07                                      \\ \hline
HGL Depth                     & 0.05            & 0.00              & 0.05                                      \\ \hline
Gas Concentrations            & 0.02            & 0.08              & 0.08                                      \\ \hline
Smoke Concentration           & 0.14            & 0.13              & 0.19                                      \\ \hline
Pressure, Closed Compartment  & 0.01            & 0.21              & 0.21                                      \\ \hline
Pressure, Open Compartment    & 0.01            & 0.15              & 0.15                                      \\ \hline
Velocity                      & 0.07            & 0.03              & 0.08                                      \\ \hline
Heat Flux                     & 0.05            & 0.10              & 0.11                                      \\ \hline
No.~Activated Sprinklers      & 0.00            & 0.15              & 0.15                                      \\ \hline
Sprinkler Activation Time     & 0.00            & 0.06              & 0.06                                      \\ \hline
Cable Failure Time            & 0.00            & 0.12              & 0.12                                      \\ \hline
Smoke Alarm Activation Time   & 0.00            & 0.34              & 0.34                                      \\ \hline
\end{tabular}
\end{center}
\label{Uncertainty}
\end{table}



\section{Calculating Model Uncertainty}

This section describes a method for calculating the {\em model uncertainty}~\cite{McGrattan:Metrologia}. Specifically, this entails developing formulae for the mean and standard deviation of a statistical distribution like the one shown in Fig.~\ref{bell_curve}. These formulae are functions solely of the model predictions and the experimental measurements against which the model is compared. The objective is to characterize the performance of the model in predicting a given quantity of interest (e.g., the hot gas layer temperature) with two parameters; one that expresses the tendency for the model to under or over-predict the true value of the quantity and one that expresses the degree of scatter about the true value.

The predicted and measured values of the quantity of interest are obtained from one or more validation studies. Figure~\ref{temp_history} is a typical example of a comparison of model and measurement. Given that usually dozens of such measurements are made during each experiment, and potentially dozens of experiments are conducted as part of a test series, hundreds of such plots can be produced for any given quantity of interest.
\begin{figure}[ht]
\begin{center}
\includegraphics[height=2.5in]{FIGURES/sample_time_history}
\end{center}
\caption[Sample time history plots]{Example of a typical time history comparison of model prediction and experimental measurement.}
\label{temp_history}
\end{figure}
Usually, the data is condensed into a more tractable form by way of a single metric with which to compare the two curves, like the ones shown in Fig.~\ref{temp_history}. Peacock et al.~\cite{Peacock:FSJ1999} discuss various possible metrics. The metric most often used in the FDS validation process is the difference between the peak values of the time-averaged predicted and measured time histories. Thus, for each measurement point, $i$, there is a pair of \underline{E}xperimental and \underline{M}odeled values, $(E_i,M_i)$, where $i$ ranges from 1 to $n$ and both $E_i$ and $M_i$ are positive numbers expressing the greatest rise in the value of the measured/predicted quantity above its ambient. Here, $n$ is the total number of measurement points of a particular quantity of interest from all the experiments.

As mentioned above, measurements from full-scale fire experiments often lack uncertainty estimates. In cases where the uncertainty is reported, it is usually expressed as either a standard deviation or confidence interval about the measured value. In other words, there is rarely a reported systematic bias in the measurement because if a bias can be quantified, the reported values are adjusted accordingly. For this reason, assume that a given experimental measurement, $E_i$, is normally distributed about the ``true'' value, $\theta_i$, and there is no systematic bias:
\be
   E \; | \; \theta \sim N(\theta \; , \; \sigma_E^2) \label{expunc}
\ee
The notation\footnote{Note that the subscript, $i$, has been dropped merely to reduce the notational clutter.} $E \; | \; \theta$ means that $E$ is conditional on a particular value of $\theta$. This is the usual way of defining a likelihood function. It is convenient to use the so-called delta method\footnote{Given the random variable $X \sim N(\mu,\sigma^2)$, the delta method~\cite{Oehlert:1992} provides a way to estimate the distribution of a function of $X$: \[g(X) \sim N \left( g(\mu) + g''(\mu) \, \sigma^2/2 \, , \, (g'(\mu) \, \sigma)^2\right)\]} to obtain the approximate distribution
\be
   \ln E \; | \; \theta \sim N \left( \ln \theta - \frac{\widetilde{\sigma}_E^2}{2} \, , \,\widetilde{\sigma}_E^2 \right) \label{eeq}
\ee
The purpose of applying the natural log to the random variable is so that its variance can be expressed in terms of the relative uncertainty, $\widetilde{\sigma}_E=\sigma_E/\theta$. This is the way that experimental uncertainties are reported. In addition, the results of past validation exercises, when plotted as shown in Fig.~\ref{scatterplot}, form a wedge-shaped pattern that suggests that the difference between predicted and measured values is roughly proportional to the magnitude of the measured value.

It cannot be assumed, as in the case of the experimental measurements, that the model predictions have no systematic bias. Instead,
it is assumed that the model predictions are normally distributed about the true values
multiplied by a bias factor, $\delta$:
\be
   M \; | \; \theta \sim N \left(\delta \, \theta \, , \, \sigma_M^2 \right) \label{mdist}
\ee
The standard deviation, $\sigma_M$, and the bias factor, $\delta$, represent the model uncertainty.
Again, the delta method renders a distribution for $\ln M$ whose parameters can be expressed in terms of a
relative standard deviation:
\be
   \ln M \; | \; \theta \sim N \left(\ln \delta +\ln \theta - \frac{\widetilde{\sigma}_M^2}{2} \; , \;
   \widetilde{\sigma}_M^2 \right) \quad ; \quad \widetilde{\sigma}_M=\frac{\sigma_M}{\delta \, \theta} \label{meq}
\ee
Combining Eq.~(\ref{eeq}) with Eq.~(\ref{meq}) yields:
\be
   \ln M  - \ln E = \ln(M/E) \sim N \left( \ln \delta - \frac{\widetilde{\sigma}_M^2}{2}+\frac{\widetilde{\sigma}_E^2}{2} \; ,
   \; \widetilde{\sigma}_M^2+\widetilde{\sigma}_E^2 \right) \label{lnMeq}
\ee
To estimate the mean and standard deviation of the distribution\footnote{The assumption that $\ln(M/E)$ is normally distributed has been tested for each quantity of interest discussed in the chapters ahead. The results are shown in Section~\ref{normality_tests}.}, first define:
\be
   \overline{\ln (M/E)} = \frac{1}{n} \, \sum_{i=1}^n \, \ln (M_i/E_i)
\ee
The least squares estimate of the standard deviation of the combined distribution is defined as:
\be
   \widetilde{\sigma}_M^2 + \widetilde{\sigma}_E^2 \approx \frac{1}{n-1} \sum_{i=1}^n \,
   \left[ \ln (M_i/E_i) - \overline{\ln (M/E)}  \right]^2 \label{stdev}
\ee
Recall that $\widetilde{\sigma}_E$ is known and the expression on the right can be evaluated using the pairs of measured and
predicted values. Equation~(\ref{stdev}) imposes a constraint on the value of the experimental uncertainty, $\widetilde{\sigma}_E$. A further constraint is that $\widetilde{\sigma}_M$ cannot be less than $\widetilde{\sigma}_E$ because it is not possible to demonstrate that the model is more accurate than the measurements against which it is compared. Combining the two constraints leads to:
\begin{equation}
   \widetilde{\sigma}_E^2 < \frac{1}{2} \mathrm{Var}\Big( \ln (M/E) \Big)
\end{equation}
An estimate of $\delta$ can be found using the mean of the distribution:
\be
   \delta \approx \exp \left( \overline{\ln (M/E)} + \frac{ \widetilde{\sigma}_M^2}{2}-\frac{\widetilde{\sigma}_E^2}{2} \right) \label{delta}
\ee
Taking the assumed normal distribution of the model prediction, $M$, in Eq.~(\ref{mdist}) and using
a Bayesian argument\footnote{The form of Bayes theorem used here states that the posterior distribution is the product of
the prior distribution and the likelihood function, normalized by their integral:
$f(\theta|M)= p(\theta) \, f(M|\theta)/\int p(\theta) \, f(M|\theta) \, d\theta$.
A constant prior is also known as a Jeffreys prior~\cite{Gelman:Stats}.}
with a non-informative prior for $\theta$, the posterior distribution can be expressed:
\be
   \delta \, \theta \; | \; M \sim N \left( M \; , \; \sigma_M^2 \right) \label{thetaeq}
\ee
The assumption of a non-informative prior implies that there is not sufficient information about the
prior distribution (i.e., the true value) of
$\theta$ to assume anything other than a uniform\footnote{A uniform distribution means that for any two equally sized intervals of the real line,
there is an equal likelihood that the random variable takes a value in one of them.} distribution.
This is equivalent to saying that the modeler has not biased the model input parameters to compensate for a known
bias in the model output. For example, if a particular model has been shown to over-predict compartment temperature, and the modeler has reduced the specified heat release
rate to better estimate the true temperature, then it can no longer be assumed that the prior distribution of the true temperature is uniform.
Still another way to look at this is by analogy to target shooting. Suppose a particular rifle
has a manufacturers defect such that, on average, it shoots 10~cm to the left of the target. It must be assumed that any given shot by a marksman without this knowledge is
going to strike 10~cm to the left of the intended target. However, if the marksman knows of the defect, he or she will probably aim 10~cm to the right of the
intended target to compensate for the defect. If that is the case, it can no longer be assumed that the intended target was 10~cm to the right of the bullet hole.

The final step in the derivation is to rewrite Eq.~(\ref{thetaeq}) as:
\be
   \theta \; | \; M \sim N \left( \frac{M}{\delta} \; , \; \widetilde{\sigma}_M^2 \left( \frac{M}{\delta} \right)^2 \right) \label{truth}
\ee
This formula has been obtained\footnote{Note that if $X \sim N(\mu,\sigma^2)$, then
$cX \sim N ( c \mu , (c \sigma)^2)$.} by dividing by the bias factor, $\delta$, in Eq.~(\ref{thetaeq}). To summarize, given a model prediction, $M$,
of a particular quantity of interest (e.g., a cable temperature), the true (but unknown) value of this quantity is normally distributed. The mean value
and variance of this normal distribution are based solely on comparisons of model predictions with past experiments that are similar to the particular fire
scenario being analyzed. The performance of the model is quantified by the estimators of the parameters, $\delta$ and $\widetilde{\sigma}_M$, which
have been corrected to account for uncertainties associated with the experimental measurements.

When computing the relative error between measured and predicted values, an additional step is performed to ensure that the bias factor, $\delta$, is not skewed by a large number of data points at any particular region in the scatter plot. The approach used for this procedure is called a regressogram, i.e., a bin-smoothed estimator function~\cite{Tukey:1}. This approach accounts for cases in which small measured values are compared to small predicted values, which can result in a large relative error. In these cases, the calculated bias factor might not be representative of the overall model bias, especially for larger measured and predicted values. Alternatively, a regressogram treats the average values throughout the scatter plot equally by subdividing the scatter plot into bins and normalizing each bin by the number of local data points. The regressogram estimator function is implemented as follows. For each scatter plot, the $x$-axis is subdivided into 10 equally spaced bins. Each bin is assigned a weight that is inversely proportional to the number of points in the bin; a bin with more points is assigned a smaller weight, and a bin with fewer points is assigned a larger weight. Finally, when the relative error is calculated, each bin is multiplied by its respective bin weight.



\section{Example}

This section describes how to make use of Eq.~(\ref{truth}). Referring to the sample problem given above, suppose a fire model is being used to estimate the likelihood that
electrical control cables could be damaged due to a fire in a compartment.
Damage is assumed to occur when the surface temperature of any cable reaches 200~$^\circ$C.
What is the likelihood that the cables would be damaged if the
model predicts that the maximum surface temperature of the cables is 175~$^\circ$C. Assuming that the input parameters are not in question, the following
procedure is suggested:
\begin{enumerate}
\item Assemble a collection of model predictions, $M_i$, and experimental measurements, $E_i$, from
past experiments involving objects with similar thermal characteristics as the cables in question.
How ``similar'' the experiment is to the hypothetical scenario under study can be quantified by way of
various parameters, like the thermal inertia of the object, the size of the fire, the size of the compartment, and so on. Obtain estimates of the
experimental uncertainty from those who conducted the experiments or follow the procedure outlined by Hamins~\cite{NUREG_1824_Sup_1}. Express the
experimental uncertainty in relative terms, $\widetilde{\sigma}_E$.
\item Calculate the bias factor, $\delta$, and relative standard deviation, $\tilde{\sigma}_M$, from Eqs.~(\ref{delta}) and (\ref{stdev}), respectively.
\end{enumerate}
Consider the distribution, Eq.~(\ref{truth}), of the ``true'' temperature, $\theta$, shown graphically in Fig.~\ref{bell_curve}.
The vertical lines indicate the ``critical'' temperature at which damage is assumed to occur ($T_c=200$~$^\circ$C), and the temperature predicted by the
model (175~$^\circ$C). Given an ambient temperature of 20~$^\circ$C, the predicted temperature rise, $M$, is 155~$^\circ$C.
The mean and standard deviation in Eq.~(\ref{truth}) are calculated:
\be \mu = 20 + \frac{M}{\delta} = 20 + \frac{155}{1.13} = 157 \; ^\circ \hbox{C}  \quad ; \quad
   \sigma = \widetilde{\sigma}_M \, \frac{M}{\delta} = 0.20 \times \frac{155}{1.13} = 27 \; ^\circ \hbox{C}  \ee
respectively. The shaded area beneath the bell curve is the probability that the ``true'' temperature can exceed the
critical value, $T_c=200$~$^\circ$C, which can be expressed via the {\em complimentary error function}:
\be P(T>T_c) = \frac{1}{2} \hbox{erfc} \left( \frac{T_c - \mu}{\sigma \sqrt{2}} \right) = \frac{1}{2} \hbox{erfc} \left( \frac{200 - 157}{27 \sqrt{2}} \right) \approx 0.06  \ee
This means that there is a 6~\% chance that the cables could become damaged, assuming that the model's input parameters are not
subject to uncertainty.


\section{Additional Considerations}


Keep in mind that for any fire experiment, FDS might predict a particular quantity accurately (within the experimental uncertainty bounds, for
example), but another quantity less accurately. For example, in the a series of 15 full-scale fire experiments conducted at NIST in 2003, sponsored
by the U.S.~Nuclear Regulatory Commission, the average hot gas layer (HGL) temperature predictions were nearly within the accuracy of the measurements
themselves, yet the smoke concentration predictions differed from the measurements by as much as a factor of 3. Why? Consider the following issues
associated with various types of measurements:
\begin{itemize}
\item Is the measurement taken at a single point, or averaged over many points? In the example above, the HGL temperature is an average of many
thermocouple measurements, whereas the smoke concentration is based on the extinction of laser light over a short length span. Model error tends to
be reduced by the averaging process, plus most fire models, including FDS, are based on global mass and energy conservation laws that are expressed
as spatial averages.
\item Is the measured quantity time-averaged or instantaneous? For example, a surface temperature prediction is less prone to error in comparison to a
heat flux prediction because the former is, in some sense, a time-integral of the latter.
\item In the case of a point measurement, how close to the fire is it? The terms ``near-field'' and ``far-field'' are used throughout this Guide to describe
the relative distance from the fire. In general, predictions of near-field phenomena are more prone to error than far-field. There are exceptions,
however. For example, a prediction of the temperature directly within the flaming region may be more accurate than that made just a fire diameter
away because of the fact that temperatures tend to stabilize at about 1000~$^\circ$C within the fire itself, but then rapidly decrease away from the
flames. Less accurate predictions typically occur in regions of steep gradients (rapid changes, both in space and time).
\end{itemize}


